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ABSTRACT 

The luminous blue variable r\ Carinae is surrounded by a complex and highly structured nebula of 
ejected material. The best-studied and axisymmetric components of this outflow consist of bipolar 
lobes (the "homunculus" ) and an equatorial "skirt." Recent proper motion measurements suggest 
that the skirt was ejected at the same time as the lobes, contrary to the assumptions of all current 
theoretical models for the formation of the nebula (which use the skirt to collimate stellar winds into 
lobes). We present a magnetohydrodynamic (MHD) stellar wind model that produces an outflowing 
disk and bipolar lobes in a single, steady-state wind. The basic model consists of a wind from a 
rotating star with a rotation-axis-aligned dipole magnetic field. The azimuthal component of the 
magnetic field, generated by stellar rotation, compresses the wind toward the equator and also toward 
the rotation axis, simultaneously producing an outflowing disk and jet. We use numerical MHD 
simulations to study the wind for various amounts of stellar rotation and to show a range of wind 
morphologies. In order to produce wide angle lobes similar to the homunculus (which have roughly 
a 30 degree opening angle), a high-speed polar wind (with enhanced energy density) from the star 
is also required. In that case, the structure of the wind bears a remarkable resemblance to the skirt 
plus homunculus morphology of the r\ Car nebulae, and a significant fraction of the stellar angular 
momentum is carried away by the wind. Although the model assumes a steady-state wind (rather 
than an eruption) and thermal wind driving (rather than radiation pressure), the structure of the 
wind is encouraging. 

Subject headings: circumstellar matter — magnetic fields — MHD — stars: individual (77 Carinae) — 
stars: winds, outflows 



1. INTRODUCTION 

During the "great eruption," beginning in 1837 and 
lasting until about 1858, the luminous blue variable 
77 Carinae ejected more than one solar mass of mate - 
rial (for a review, see iDavidson fc Humphreys! Il9971. 
That ejected mater ial has been well studied (e.g., 
IDavidson et al.l2001[) and is highly structured, consisting 
of an outflowing equatorial "skirt" and bipolar lobes (the 
hourglass-shaped "homunculus" ) . The presence of the 
bipolar lobes has been exp l ained b y some authors (e.g., 
i Frank. Balick. fc Davidsonl 119951 iDwarkadas fc Balick 
1998; Laneer. Garcfa-Segura, & Mac Low 1999) with in- 
teracting stellar wind (ISW) models. In such models, 
a previously ejected and slowly moving equatorial flow 
(presumably related to the skirt) acts as an inertial bar- 
rier that channels a subsequent flow to ward the poles, 
formi ng the homunculus (but see also iGonzalez et al.l 
2004). 

However, recent proper motion measurements by 
iMorse et alJ (|2P01) reveal that much of the material in 
the skirt has the same dynamical age as the lobes. These 
observations suggest that the skirt was ejected at the 
same time as the lobes (within the uncertainty of a few 
years) in 1846, during the great eruption, and thus rule 
out the ISW models. Clearly, there is a need for a theo- 
retical model that produces bipolar lobes and an equato- 
rial flow, simultaneously, in a single wind. As far as we 
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know, no such model yet exists in the literature. 

An alternative to the ISW scenario is shaping of the 
outflow by some sort of magnetohydrodynamic (MHD) 
proc ess. This ex p lanat ion has observational support 
from lAitken et all l)1995|) . who measured the polariza- 
tion of mid-infrared emission from the spatially-resolved 
homunculus. Their explanation for the structure and 
degree of polarization requires a magnetic field within 
the homunculus with a strength in the range of milli- 
Gauss. The homunculus is surrounded by material that 
was ejected prior to the great eruption, so the magnetic 
field in the homunculus cannot be the result of swept-up 
and compressed ISM field, and it must have originate d 
from the central source of the wind ijAitken et alJ ll995). 

The influence of magnetic fields on stella r winds has 
been studied by many authors. Sakurai ( 1983) ex ~ 
tended the magnetohy d rodyn amic (MHD) wind the- 
ory of iWeber fc Davis! l)1967|) and studied isotropic, 
thermally-driven stellar winds, embedded in a radial 
magnetic field. Sakurai showed that when the star ro- 
tates, the magnetic field is twisted azimuthally, which 
results in a collimation of the wind toward the r otation 
axis, as well as an increased wind s peed (see alsolMichel 
1969: IBlandford fc Pa.vnei 11982^. IWa,shimi fc ShiCT 
(|1993i hereafter WS93) considered the effect of a dipole 
magnetic field on an otherwise isotropic, pressure-driven 
stellar wind. They showed that, when the star ro- 
tates rapidly, the azimuthally twisted dipole produces 
the same collimation shown by Sakurai, but also results 
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in an enhancement of the wind along the magnetic equa- 
tor. In other words, the wind from a rapidly rotating 
star with an aligned dipolar magnetic field contains both 
an outflowing disk and a jet. 

It seems plausible that the skirt in the r\ Car nebula 
was produced by a similar mechanism, as the skirt may 
correspond to the outflowing disk in such a wind. On the 
other hand, the homunculus consists of wide angle lobes, 
and does not match the narrow jet morphology predicted 
by the standard, rotating MHD wind theory. So the idea 
of an isotropic wind driving mechanism, plus a rotating 
dipole magnetic field, must be modified to produce a 
wider angle polar flow. In this work, we follow the work 
of WS93 and explore the possibility that the skirt and 
homunculus were formed during the great eruption in a 
single wind from a rotating star with a dipole magnetic 
field. In order to produce the wide-angle, bipolar lobes of 
77 Car, we find that the wind driving must have been more 
energetic near the poles than at lower latitudes. Section 
[21 contains a description of the general model and the 
assumptions therein. 

Using numerical, 2.5-dimensional (axisymmetric) , 
MHD simulations, we follow the acceleration of the stel- 
lar wind from the surface of the star to beyond 100 stellar 
radii, self-consistently capturing the interaction between 
the stellar magnetic field and the wind plasma. In order 
to understand each physical process active in the wind 
and to demonstrate a range of wind morphologies, we 
carry out a limited parameter study. Section contains 
a description of the simulations and the parameters we 
explored. Section0|contains the results of the parameter 
study, which are quite general, and may be applicable 
to a variety of structured stellar outflows, particularly 
those with quadrupolar symmetry, such as pulsar winds 
and proto-planetary nebulae. One of the simulations pro- 
duces in a steady-state wind with structures that bear 
a remarkable resemblance to the skirt and homunculus 
morphology of r\ Car. 

2. GENERAL MODEL AND PHYSICAL PARAMETERS 

Axisymmetric outflow nebulae are the result of a coher- 
ent but poorly understood process that probably arises 
in very close proximity to the stellar surface. Thus, one 
of the ultimate goals of any wind study is to link observed 
wind parameters (and morphology) to the conditions at 
or very near the source of the wind. In this paper, we 
explore the possibility that the outflow from 77 Car dur- 
ing the great eruption was produced by a wind from an 
isolated star. We hope to constrain the possible condi- 
tions near the surface of r\ Car during the great eruption. 
Knowing these conditions, and how they change in time, 
is part of the bigger picture of understanding the internal 
structure and evolution of extremely massive stars. 

We shall assume that 77 Car had a dynamically impor- 
tant magnetic field on its surface during outburst. Mag- 
netic fields are capable of producing structure in other- 
wise isotropic flows, particularly in systems with rota- 
tional energy, where the twisting of magnetic field lines 
naturally leads to a collimation of the flow along the ro- 
tation axis. We will employ an MHD formulation of the 
problem, adopting all of the assumptions therein. This 
assump tion is suppo r ted, i n principle, by the observa- 
tions of lAitken etaD l)1995|) . 

In general, any MHD stellar wind is characterized by 



several key components: magnetic field geometry, mag- 
netic field strength, stellar rotation rate and gravity, and 
the wind driving mechanism. These components may 
have time dependence as well as spatial variations along 
the surface of the star. In addition, they are all likely 
to be interrelated. For example, the rotation rate may 
influence a stellar dynamo, which determines the field 
strength and geometry, which in turn, may influence the 
structure and temperature in a corona, etc. However, 
since there is currently no comprehensive theory relating 
all key components with one another, we treat them as 
independent parameters. 

Following the work of WS93, we consider pressure 
(thermal) driving for the wind, using a polytropic (7 = 
1.05) equation of state, as appropriate for solar wind 
studies. The assumption of thermal driving is not very 
appropriate for 77 Car, as radiation pressure is likely 
to h ave been more important du ring the great eruption 
('e.g.. lDwarkadas fc Owo cki 2002). However, since this is 
an initial study, our assumption is for simplicity, and the 
general results will be valid for any wind driving mech- 
anism. Also as in WS93, the magnetic field anchored in 
the surface of the star is a rotation axis aligned dipole 
with a fixed field strength. A dipolar magnetic field 
is the most reasonable choice, as it is the lowest order 
multipole occurring in nature and reflects the symme- 
tries inherent to a rotating star. Even the sun's mag- 
netic field is dominated by the dipole component out- 
side a few solar r adii over most of the magnetic cycle 
l|Bravo et al.lll998lK Significant progress has been made 
in recent years toward understanding the effe ct of dipole 
fields on stellar win ds (e.g., WS9 3; Babel & Mon tmerlg 
19971 iPorterl Il997t iKeppens fc Coedbloedl 119991 1200(1 
Matt et all200rtlud-no 1 ila fc Owockl2002D. and this lit- 
erature provides a background for the work contained 
here. 

There are four fundamental parameters, one for each 
of the energies in the system, which can be written in 
terms of characteristic speeds. The thermal energy is 
a measure of the wind driving (via pressure gradients) 
and is represented by the sound speed at the surface of 
the star, c s . The gravitational potential energy must 
be overcome by the wind and is parameterized by the 
escape speed from the surface, v csc . The rotational ki- 
netic energy is represented by the rotation speed at the 
equator, w ro t, and the magnetic energy by the Alfven 
speed at the equator, wa- The conditions in the wind, 
including the location of various critical points (where 
the wind flow speed equals other characteristic speeds), 
outflow speed, mass flux, etc., will be completely deter- 
mined, self-consistently, from the conditions present on 
the stellar surface (i.e., at the base of the wind). 

For simplicity, and so that we can first understand the 
shaping of individual winds, we only consider steady- 
state wind solutions. Below, we present winds with fea- 
tures that resemble both the skirt and homunculus, but 
since the 77 Car nebula was ejected in an outburst, we 
are only tackling part of the problem. Specifically, we as- 
sume that our steady-state solution represents the wind 
from 77 Car during the outburst. Since the true outburst 
lasted for a few years (or decades, at most), the neb- 
ula's evolution from outburst until the present day was 
likely to have been influenced by previous and subsequent 
winds. Thus, our steady-state models do not address 
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the conditions responsible for or resulting from the time- 
dependent outburst phenomenon as a whole. However, 
since there are currently no other stellar wind models 
that produce a disk and lobes simultaneously, we hope 
that the present work can serve as a starting point for 
more complex models, and that it may be more widely 
applicable to other classes of structured stellar outflows. 

The rotation rate of 77 Car is unknown, so we con- 
sider a full range of rotation rates in our models, from 
no rotation to a significant fraction of breakup speed. 
Similarly, there are currently no direct measurements 
of a magnetic field on the surface of 77 Car. Further- 
more, the possible field strength on the star during the 
great eruption, which is relevant here, is completely un- 
constrained by observations. However, we assume that 
the field was strong enough to dynamically influence the 
wind, and we pick a single value for the field strength. 
We choose not to vary the magnetic field, for simplic- 
ity, and since the eff ect of various dipole field stre ngths 
has been studied by IKe ppens fc Goedblo'edl <|200Cf) and 
lud-Doula fc Owockil |2002). In general, for the magnetic 
field to be dynamically important, the magnetic energy 
density must be comparable to the kinetic energy den- 
sity in the wind. This energy balance gives a rough lower 
limit on the magnetic field strength on the stellar surface 
of 
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where we have assumed that, during the outburst, 77 
Car had a radius of i?* = lOOi? and drove a wind 
with a mass outflow rate of M w — 1O _1 M yr _1 and 
a speed of tJqq = 600 km s -1 . These assumptions 
are consistent with observatio nally determined values 
ijDavidson fc Humphreys! 11997(1 . For the solar wind, 
equation^gives a minimum field of 0.1 Gauss, though the 
sun maintains a 1-2 Gauss dipole magnetic fie ld through- 
out most of the solar cycle l|Bravo et al.l |l998) . By draw- 
ing upon the example of the sun, we choose a dipole field 
strength of 2.5 x 10 4 Gauss, measured at the equator of 
77 Car. This is a strong field for a giant star, but it is re- 
quired if magnetic shaping is to be important. Note that 
the wind from 77 Car during outburst was extremely pow- 
erful, compared to other classes of stellar winds. If the 
process that produces an increase in wind energy during 
outburst also produces an increase in dynamo activity, 
it may be natural to assume that the magnetic field will 
also exist in a heightened state during outburst. 

We know from the work of WS93 that the effect of a 
rotating dipolar field on an otherwise isotropic wind is 
to, simultaneously, compress the flow toward the mag- 
netic equator and collimate it toward the rotation axis. 
The resulting outflowing disk superficially resembles the 
morphology of the skirt around 77 Car, but the colli- 
mated jet is too narrow to explain the wide-angle, bipolar 
lobes that comprise the homunculus. In order to produce 
these lobes, the wind driving could be more energetic 
near the poles than at lower latitudes — a similar modi- 
fication to the rota ti ng di pole wind model was used by 
iTanaka fe^ Washimil l|2002|) to reproduce the three-ring 



structure of SN 1987A. This enhancement of the po- 
lar wind effectively reduces the dynamical influence of 
the magnetic field at high latitudes, reducing the col- 
limation in that region and resulting in a polar flow 
with a wide opening angle. To this end, we present 
some models with an increased temperature on the stel- 
lar poles, and since we are considering a thermally driven 
wind, the higher temperature results in a more energetic 
flow from the poles. There is justification for an en- 
hanced polar wind from 77 Car. First, radiatively driven 
winds from rapidly-rotating, luminous blue variable stars 
are e xpected to be more energe tic in the polar direc- 
tion llDwarkadas fc Ow ockil 200211. Second, ob servations 
of rj Car ijSmith et all 120031 Ivan Boekel et al J 12003]) re- 
veal that the present-day stellar wind is significantly en- 
hanced in the polar direction. Furthermore, for more 
general support, observations of the solar wind show that 
the kinetic energy and outflow speed is larger (by a factor 
of ~ 2 ) above 30 degrees la titude than at lower latitudes 
(e.g.. lGoldstein et al.lll996jh 

In summary, we consider steady-state, MHD winds 
from a star with a rotation-axis-aligned dipole magnetic 
field and varying amounts of rotation. The wind is ther- 
mally driven, and we consider two possibilities: cases 
with isotropic wind driving, and cases with an enhanced 
polar wind. The resulting winds (described in J|J) display 
an encouraging range of morphologies. The case with the 
fastest rotation and with an enhanced polar wind most 
resembles the 77 Car nebula. 

3. MHD SIMULATION DETAILS 

Here we present numerical MHD simulations of sev- 
eral models, each with a different set of stellar surface 
values of c s , tj csc , u ro t, and i>a- It is not our intention 
to carry out a complete parameter study, which should 
include variations of all of these parameters, plus vari- 
ations of the magnetic field geometry, and wind driving 
anisotropy. Instead, our aim is to determine what sur- 
face conditions may be required to reproduce the sort of 
structures seen in the 77 Car outburst (i.e., the skirt and 
homunculus). 

The cases we ran were chosen to illustrate, in a step- 
by-step manner, the individual effect of each physical 
process, leading up to a "final" model that resembles the 
77 Car nebula. In each case, the parameters are chosen 
and held fixed on the stellar surface, and the simula- 
tions run until the entire computational domain is filled 
with a steady wind emanating from the stellar boundary. 
This method gives the steady-state solution of the flow 
within the computational grid and determined solely by 
the conditions held fixed on the stellar surface boundary. 

3.1. Numerical Method 

We use the 2.5-dimensional MHD code of iMatt et alJ 
(2002), which we describe briefly here. The reader will 
find fur ther details of the code in IMatt et al.l J2002L an d 
also in iGoodson. Winglee. fc Bohml 119971: lMatT f2002). 
which solves the ideal (non-resistive) MHD equations 
using a two-step Lax-Wend roff, finite difference scheme 
(Richt mver fc Mortonl lT967) in cylindrical (r, 0, z) geom- 
etry. The formulation of the equations allows for a poly- 
tropic equation of state (we use 7 = 1.05), includes a 
source term in the momentum and energy equations for 
point-source gravity, and makes the explicit assumption 
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of axisymmetry {d /dtp = for all quantities). The fun- 
damental plasma quantities in the equations are the vec- 
tor magnetic field B, mass density p, vector momentum 
density pv (where v is the velocity), and energy density 
e = O.bpv 2 + P/(7 — 1), where P is the thermal pressure. 
Due to the geometry and symmetry in the system, it is 
often useful to decompose the vectors into the azimuthal 
and poloidal components. The azimuthal component is 
that in the cylindrical <f> direction and denoted by a sub- 
script "</>," and the poloidal component is that in the r-z 
plane denoted by a subscript "p" (e.g., v = v p f> + v^cf) 
and v p f> = v r f + v z z). 

The computational domain consists of four nested, 
grids (or "boxes") in the cylindrical r-z plane. Each box 
contains 401 x 400 gridpoints (in r and z, respectively) 
with constant grid spacing. The boxes are nested concen- 
trically, so that the inner box represents the smallest do- 
main at the highest resolution. The next outer box rep- 
resents twice the domain size with half the spatial resolu- 
tion (an so on for other boxes). With four boxes, the to- 
tal computational domain is eight times larger than that 
of the innermost grid, but requiring only four times the 
computational expense. This computational efficiency 
allowed us to run the several cases necessary for the cur- 
rent study. A circular boundary, centered on the origin 
and with a radius of 30 gridpoints, represents the surface 
of the star. Assuming R* — 100i?© for ij Car during out- 
burst, the innermost grid then has a spatial resolution of 
roughly 3.33i?0, and a domain size of 13. 3i?*. Similarly, 
the outermost (fourth) box has a resolution and domain 
size of 26.7-R© and 107-R*, respectively. 

3.2. Boundary and Initial Conditions 

We use standard outflow conditions on the box bound- 
aries, appropriate for wind studies. Namely, along the 
outermost boundary in r and z in the fourth box, a con- 
tinuous boundary condition (in which the spatial deriva- 
tive across the boundary is zero for all quantities) allows 
the stellar wind to flow out of the computational domain 
undisturbed. In all grids, we enforce reflection symmetry 
across the equatorial [z = 0) plane. Along the rotation 
axis (r = 0), we require the azimuthal and radial com- 
ponents of magnetic field and velocity to be zero, and all 
other quantities are equal to their value at r = dr (where 
dr is the grid spacing). 

The surface of the star is represented by a spherical 
inner boundary, centered at the origin (r = z = 0). This 
boundary is the source of the wind that is to fill the 
computational domain, and so the conditions here en- 
tirely determine the solution of the system. Following the 
wind from the very surface of a star and self-consistently 
capturing the interaction between the wind plasma and 
a magnetic field that is anchored into the rotating stel- 
lar surface is a complex problem. In a real star, the 
properties of the wind (mass flux, velocity, etc.) are de- 
termined by the balance of forces resulting from (e.g.) 
thermal pressure gradients, gravity, centrifugal forces, 
and through interaction with the poloidal magnetic field 
(which distorts when pushed and also pushes back) and 
the azimuthal magnetic field — which is generated by a 
twisting of the poloidal field as outflowing plasma tries 
to conserve its angular momentum. 

In order to capture these interactions within the frame- 
work of our second-order finite difference scheme, we em- 



ploy a four-layer boundary for the star, on which the var- 
ious parameters are set as follows: For all gridpoints such 
that R < 34.5 (where R = (r 2 + z 2 ) 1 / 2 in units of the 
grid spacing in the innermost box) , the poloidal velocity 
is forced to be parallel with the poloidal magnetic field 
(v p || B p ). Where R < 33.5, p and P are held constant 
(in time) at their initial values. For R < 32.5, v p is held 
at zero, while v^, is held at corotation with the star. For 
R < 31.5, Bp field is held at its initial, dipolar value, 
while B^ is set so that there is no poloidal electric cur- 
rent at that layer (which gives it a dependence on the 
conditions in the next outer layer, 31.5 < R < 32.5). We 
consider the spherical location R = 30 to be the surface 
of the star, as this is where all quantities are held fixed 
and is thus the absolute base of the wind. 

These boundary conditions capture the behavior of a 
wind accelerated thermally from the surface of a rotat- 
ing magnetized star, as follows: There is a layer on the 
stellar boundary (i? > 32.5) outside of which the velocity 
not fixed, but is allowed to vary in time. In this way, the 
wind speed and direction is determined by the code in 
response to all of the forces. By holding P fixed at its 
initial value for all R < 33.5, we constrain the pressure 
gradient force (thermal driving) behind the wind to be 
constant. Also, holding the density fixed at R < 33.5 
allows the region from where the wind flows to be in- 
stantly replenished with plasma. Thus, the base of the 
wind maintains a constant temperature and density, re- 
gardless of how fast or slow the wind flows away from 
that region. The existence of a layer in which v p = 
and Bp can evolve (namely, at 31.5 < R < 32.5) allows 
Bp (and v p ) to reach a value that is self-consistently de- 
termined by the balance of magnetic and inertial forces 
(and is also necessary to maintain a negligible V • B). 
We set the poloidal velocity parallel to the poloidal mag- 
netic field for the next two outer layers (which ensures a 
smoother transition from the region of pure dipole field 
and zero velocity to a that with a perturbed field and 
outflow) for reasons of numerical stability. Setting B^ 
so that the poloidal electric current is zero inside some 
radius ensures that the field is completely anchored in a 
rotating conductor (the surface of the star). This way, 
B^ evolves appropriately outside the anchored layer ac- 
cording to the interaction with the wind plasma. 

We have extensively tested these conditions for a wide 
range of all system parameters, and have been able to 
reproduce the work of other authors (e.g., the results of 
IKeppens fc Goedbloedlll999l and WS93). The advantage 
of our boundary conditions is in their versatility — they 
capture the expected behavior for a wide range of con- 
ditions. Of particular note, they produce the correct hy- 
drodynamic solution for the case with B = 0, and give 
the same solution for a similar case with a weak (dynami- 
cally insignificant) dipole magnetic field — in which dipole 
field lines are stripped open to a radial, spli t monopole 
config uration by the spherical wind (see, e.g. JMatt et alJ 
2000). In cases with a strong dipole, closed magnetic 
loops near the equator overcome thermal and centrifugal 
forces to "shut off" the flow, resulting in a magnetically 
closed region with v p = and that corotates with the 
star (i.e., self-consistently f orming a "dead zone" ; see 
IKermeus fc C T oedbloedl l2l1o?il). 

For the simulations with enhanced polar winds (see 
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M3.3|) . the pressure at high latitudes (above 74°) on the 
star is held fixed at twice the pressure at the same spheri- 
cal radius but at lower latitudes. Formally, for all regions 
on the star such that R < 33.5 and z/r > 3.5, the pres- 
sure is doubled. The resulting factor of two pressure dis- 
continuity across 74° latitude on the star is only slightly 
larger than the typical radial pressure gradients near the 
star, in which, at the resolution of our simulations, the 
variation in pressure between gridpoints is typically a fac- 
tor of 1.3. Furthermore, the code has been tested, and 
behaves well, under shock conditions, in which disconti- 
nuities can be much larger than a factor of two. Thus, 
the enhanced polar wind does not present a difficulty for 
the numerical method. 

Ultimately, the conditions in the steady-state wind 
that emanates from the stellar surface depends only on 
the parameters held fixed on the star. However, an initial 
state for the system must be specified on the entire simu- 
lation grid, so we initialized the grid with the spherically 
symmetric wind solution of Parker (1958). As discussed 
above, the pr e ssure gradient and density given by this 
initial iParkerl l)1958f) solution is held fixed for all time 
on the star, though the velocity is allowed to vary to re- 
spond to the rotation and magnet ic field of the star. This 
appro ach is similar to the work of lKeppens &: Goedbloed| 
(1999). For the cases with an enhanced polar wind, the 
pressure is doubled on the polar region of the star, so the 
fixed pressure gradient is also double. 

3.3. Parameter Space Explored 

Our approach is to illustrate, in a logical progression 
of cases, the individual effects of each additional param- 
eter. To this end, we ran the 12 simulation models listed 
in tabled The cases named "ISO," have isotropic ther- 
mal wind driv ing, in which the pressure from the initial 
IParkerl l)1958h solution is held fixed at all latitudes on 
the star. The number in the last part of the name cor- 
responds to the rotation period in days, assuming the 
fiducial radius for 77 Car during outburst used in equa- 
tion n] The rotation period is also listed in the second 
column of tabled The cases named "HTC" are those in 
which the thermal wind driving is enhanced on a "high 
temperature cap" on the star. They are identical to the 
corresponding ISO cases, except that for all regions on 
the star above a latitude of 74°, the pressure (P) is set 
and held at double the value of the ISO case. Thus, the 
star has a polar cap with twice the temperature (and 
twice the pressure gradient and thermal energy) as for 
lower latitudes, which results in a more energetic flow 
from the polar region. 

In all models, the mass and radius of the star are as- 
sumed to be 100M Q and lOOi?©, respectively, which cor- 
responds to an escape speed of i> csc = 617 km s . The 
sound speed (= yJjP/p) on the surface of the entire 
star in the ISO models and below 74° latitude for the 
HTC models is c s = 137 km s _1 (and so c s = 194 km 
s _1 in the polar region of the HTC models). The values 
of c s and tj csc determine the solution to the velocity at 
all radii in a (non-rotating, non-magnetic) hydrodynamic 
wind. We chose this value of c s so that the outflow speed 
reaches several hundred km s _1 , far from the star, and 
we set our density normalization (which is the same for 
all models) such that the mass outflow rate is around 
10 _1 AfQ yr -1 , appropriate for r\ Car during outburst. 



TABLE 1 

Simulation Parameters' 1 



Case 


T 


Vrot C 
»asc 


"A c 

^GSC 




(days) 






ISO°oo 


00 








ISO00 


00 





0.35 


ISO2000 


2000 


0.004 


0.35 


ISO200 


200 


0.041 


0.35 


ISO100 


100 


0.082 


0.35 


ISO50 


50 


0.165 


0.35 


HTC°oo 


00 








HTCoo 


00 





0.35 


HTC2000 
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a The value of c B /v CBC is the same in all models and equal to 
0.223; see text. 

b The "ISO" models have isotropic wind driving, while the 
"HTC" models include a high temperature polar region on the 
star. 

c Given is the value at the equator and surface of the star. 



The magnetic field, when nonzero, is always a rotation- 
axis-aligned dipole with a strength of i?» = 2.5 x 10 4 
Gauss on the equator of the star (and twice that on 
the pole). With our density normalization, this gives an 
Alfven speed (= B/^/Anp in cgs units) on the surface of 
the star at the equator of va = 218 km s _1 . We explore 
a full range of stellar rotation up to a rotation of 23% 
of breakup speed. Thus the equatorial rotation speed on 
the star for different models ranges from w rot = to 102 
km s^ 1 . 

The ratios of v Iot /v csc and va/v csc at the equator and 
surface of the star are listed in the third and fourth 
columns of tabled for each model. The ratio of c s /v csc 
is 0.223 for all models (except on the polar region of 
the HTC models, where c s /v csc — 0.315). In order to 
establish a baseline for the ISO cases, we first ran a 
purely hydrodynamical case (B = 0) with no rotation, 
called ISO°oo. Then, to show the effect of a 2.5 x 10 4 
G dipolar field, we ran the ISOcx) case, also with no 
rotation. Next, to illustrate the twisting of magnetic 
fields, but before the twisting is dynamically important, 
we ran the ISO2000 case, in which the rotation is so slow 
that the solution is virtually identical to the ISO00 case 
(see 33}. The effects of dynamically significant rotation 
are demonstrated with the ISO200, ISO100, and ISO50 
cases. Finally, we ran an HTC model corresponding to 
each of the ISO models, so that we could see the effect 
of the anisotropic thermal wind driving for each set of 
parameters. 

The results of these simulations (presented in £|I} are 
scalable to any star with the same values of v m t/v esc , 
VA/v esc , and c s /v esc . Note that the rotation period (and 
any physical times) should then scale with R*/v, where 
v represents any of the characteristic speeds in the sys- 
tem (e.g., v esc ). In section l4~2l we calculate the outflow 

rates of mass (M w ), angular momentum (J w ), and en- 
ergy (-E w ), which can also be scaled. The mass outflow 
rate scales with vR~j. times the density normalization, 
J w oc M w vR*, and E w cx M w v 2 . Finally, the magnetic 
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field strength (B) scales with v times the square root of 

the density normalization (so that B oc \J M w v/R*). As 
an example, the sun has the same v CB c and roughly the 
same coronal sound speed and va as our models. Using 
M w ~ 10 _14 M© yr _1 , for the sun, scaling down M w and 
i?» in our models corresponds to a dipole field strengt h 
of ~ 1 Gauss, which is appropriate ( Bravo et al.l ll998h 
Thus the ISO2000 model is similar to the sun, though 
with a rotation period of 20 days. Since the solar rota- 
tion does not have a significant effect on the solar wind 
(see, e.g., WS93), we expect the ISO2000 case to dis- 
play similar global properties as the solar wind (though 
the solar wind is further modified by anisotropic coronal 
heating, high order magnetic fields, etc.). 

3.4. Accuracy of the Numerical Solution 

After the start of the simulations, a flow emanates from 
the surface of the star. In all cases, that flow accelerates 
to a speed faster than all information carrying waves (i.e., 
slow, Alfven, and fast magnetosonic waves) within the 
computational domain, usually within the first or second 
grid. To obtain a numerical steady-state solution, we ran 
each case for ~ 2 years of physical time. Since it takes 
only ~ 0.6 yr for the wind, traveling with a typical ve- 
locity of 400 km s -1 , to cross our largest simulation box, 
material that leaves the surface of the star at the begin- 
ning of the simulation has traveled well beyond the entire 
domain by the end. A measure of the fractional rate of 
change of each of the quantities p, P, v r , v^,, v z , B r , B^, 
and B z at every gridpoint quantifies the steadiness of 
the numerical solution. In the largest box, for each case, 
the average gridpoint was changing by less than 1% per 
wind crossing time by the end of the simulations. Since 
the computations require roughly 16,000 time-steps per 
crossing time, this corresponds to a fractional change of 
less than 10 -6 per time-step. As another simple test of 
the steadiness and conservation of our numerical solu- 
tions, we calculated the outflow rates of M w , J w , and 
E w (see E14.2f> at all radii within the computational do- 
main and found that these global quantities are constant 
(i.e., no variation with radius), within 3%, for all cases. 
Thus, our solutions are sufficiently steady for the present 
study. 

One general physical test of numerical MHD solutions 
is to quantify the divergence of B, which should al- 
ways and everywhere be zero. The initial state of the 
system has V • B = 0, and the code maintains this 
physical requirement to second order in space and time. 
However, numerical errors do lead to non-zero values of 
V • B, and we deal with these errors in two ways. First, 
the code solves the momentum equation using a non- 
conservative formulation of the magnetic force that pre- 
vents numerical V • B errors fro m producing unphysical 
forces l|Brackbill fc Barn es 1980). Second, we check that 
these numerical errors do not contribute significantly to 
the structure of the magnetic field, by computing the ra- 
tio of V • B to \VB\ — that is, the magnitude of the gra- 
dient of the magnitude of B — at every gridpoint. This 
ratio determines the maximum possible contribution of 
non-zero V • B to any spatial change in the magnetic 
field. For the vast majority of gridpoints, in all cases, 
this ratio is always much less than 1%. Near the stellar 
boundary, this can be larger, though never exceeding 5% 



in all magnetic cases. These V • B errors are acceptable. 

For a more comprehensive check on th e precision of our 
solutio ns, we follow the error analysis of lUstvugova et all 

(1999) , who list five integrals of motion K, A, Q, S, and 
E (their eqs. 1-5), corresponding to the conservation of 
mass, angular momentum, helicity, entropy, and energy, 
respectively. Under the conditions of ideal, axisymmet- 
ric, steady-state MHD, these integrals should be constant 
along a given magnetic field line, and thus serve as a 
stringent test of our numerical solution. In every mag- 
netic model, we calculated K, A, Q, S, and E along two 
reference field lines: one that connects to the stellar sur- 
face at 65° latitude, and one that connects to 80°. In this 
way, we sampled two latitudes in the flow, one of which 
is within the hot polar cap of the HTC models. In all 
cases, all five integral functions are conserved to within 
5%, along both reference field lines. 

For an even more stringent test, we compared the ab- 
solute value of the conserved quantity f2 to the known 
value at the stellar surface, which equals the constant 
angular ro tation rate of the star. Follo wing the error 
analysis of lKeppens fe G ocdblocd (2000), we calculated 
fi for all gridpoints in our domain, for all of the magnetic 
models. Over 90% of the domain, the error in ranges 
from 2-7%, for all cases, except one. This exception is 
the ISO50 case, in which this error is 12%. A larger error 
in f2 exists in all cases over an area of less than 10% of the 
domain that covers a narrow region near the equator with 
an opening angle of less than 10°, as seen by the star. 
This error occurs in the portion of the flow that travels 
along the lowest latitude open field lines near the dead 
zone (see @. In that portion of the flow, the maximum 
error in ranges from 20-40% in all cases, except one. 
This exception is the HTC2000 case, in which the maxi- 
mum error in this portion of the f low reaches 80%. These 
errors are comparable to those of lKeppens fe Goedbloedl 

( 2000) , who found errors in Q reaching 40% in the same 
portion of the flow as where our largest errors occur. 

4. SIMULATION RESULTS 

For the cases with nonzero magnetic field, the wind 
is significantly altered from the purely hydrodynam- 
ical solution. The behavior of the dipole magnetic 
field in a wind and it's effect on the thermally driven 
flow is discussed by man y authors (see especially 
iKeppens fc Goedbloedl 2000'). In particular, the dipole 
field lines at high latitudes are nearly radial, and thus 
wind can flow along them with little impedance. Fur- 
ther, they reach to distances far from the star, where the 
field becomes very weak. As the dipole field energy den- 
sity falls off with R~ e , and the kinetic energy density in 
the wind falls off more slowly than R~ 2 (for an acceler- 
ating wind) these high latitude field lines will open up as 
they are dominated by and carried out in the wind. On 
the other hand, at lower magnetic latitudes, the dipole 
field lines are more perpendicular to the flow, and they 
do not reach as far away from the star. If the field is 
strong enough to counteract the outward thermal pres- 
sure forces, the flow will be shut off. This produces a 
region of closed field lines around the magnetic equator, 
a "dead zone," that corotates with the star and from 
which no wind can flow (in the absence of magnetic dif- 
fusion) . 

If the star rotates, wind material tries to conserve its 
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Fig. 1. — Greyscale of (white is zero, black is strong) with selected magnetic field lines overdrawn in grey for two cases where the 
star has a constant temperature across its surface. The dashed line indicates the radial Alfven surface, and the dotted line is the sonic 
surface. The left panel contains data for the very weakly rotating case (T = 2000 days), while the right panel corresponds to faster rotation 
(T = 100 days). The dash-dotted line in the left panel is the location of the sonic surface for a simulation without rotation, and without a 
magnetic field. The light grey field line in the right panel intersects the stellar surface at 65° and indicates the path s followed in Figure 

m 



angular momentum as it travels outward so that its an- 
gular velocity decreases with distance from the rotation 
axis. When wind plasma is connected to the stellar sur- 
face by magnetic fields, the difference in angular rota- 
tion rate between the star and wind twists the field az- 
imuth ally, genera ting a component to the field (see, 
e.g.. lParkerJll958i) . If the stellar rotation is fast enough, 
magnetic forces associated with significantly influ- 
ence the flow (e.g., WS93). These effects will be further 
discussed below. 

4.1. Wind Acceleration and Shaping 

The left panel of Figure ^ illustrates the conditions in 
the steady-state wind for the ISO2000 case. Shown is a 
greyscale image of B^ (black represents 57 Gauss, and 
white is zero), the poloidal magnetic field lines (solid 
lines), and poloidal velocity vectors (arrows). The data 
is from the innermost (highest resolution) computational 
box. The star is at the lower left and has a radius of i?* . 
Poloidal magnetic field lines are drawn so that they origi- 
nate at roughly evenly spaced intervals in latitude on the 
star (their spacing is not proportional to field strength). 
We have also drawn a field line that marks the transi- 
tion between closed and open field regions (i.e., the line 
that encloses the dead zone). Note that, under steady- 
state, ideal MHD conditions (as we have here), magnetic 
field lines are also in the same direction as plasma flow 
streamlines. The rotation is so slow in the ISO2000 
case that the resulting wind is identical (within a few 
percent) to the ISO00 case, except that B^ is exactly 
zero for ISO00. Thus, the left panel of Figure ^ simply 
demonstrates the effects of a 2.5 x 10 4 G dipole magnetic 
field on an otherwise completely isotropic flow. The re- 



sults here are comparable to other similar studies (e.g., 
iKeppens fc Goedbloedl2000() . and they form the base to 
our understanding of the effects of additional processes 
in the wind (i.e., rotation and enhanced polar winds). 

In the left panel of Figure ^ the dead zone covers a 
significant fraction of the stellar surface near the equator 
and extends beyond 5 i?* above the star. Thus, all of the 
wind from the star originates from the open field region 
at high latitudes. There is flow at low latitudes, outside 
several stellar radii, but (as indicated by the field lines 
and velocity vectors) that flow originates from the polar 
region of the star. In other words, the wind from the 
poles diverges very rapidly near the star and tends more 
toward radial divergence far from the star. The flow ac- 
tually converges toward the equator within roughly 8i?» , 
as the open field lines curve around the dead zone. The 
dotted line in the left panel of Figurcn]marks the location 
of the sonic surface, where the wind velocity equals the 
local sound speed. For reference, the dash-dotted line 
marks the sonic surface for the ISO°oo case, in which 
there is no magnetic field, and the wind is completely 
isotropic. The effect of a dipole field is to shift the sonic 
surface further from the star near the equator. This is 
because the convergent flow on the equator leads to a 
shallower radial pressure gradient, and the wind is accel- 
erated more slowly there. Conversely, the sonic surface 
is moved closer to the star near the pole because the 
flow diverges faster than radial divergence, making the 
pressure gradient steeper. The dashed line marks the lo- 
cation of the poloidal Alfven surface, where the poloidal 
wind velocity equals the local poloidal Alfven speed (cal- 
culated from the B p component of the field only) . The 
Alfven surface is furthest from the star along the pole, 
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since the magnetic field is strongest there. Due to the 
opening of field lines in the wind, the poloidal magnetic 
field strength goes to zero on the equator (across which 
there is a direction reversal of the field, supported by a 
current sheet), leading to a "cusp" in the Alfven surface 
there. 

The grey scale of in the left panel of Figure ^ shows 
that the dead zone is corotating with the star, since 
B^ « there. However, in the open field region, the 
field lines are "trailing" the star as it rotates beneath 
the expanding wind, and so B^ ^ 0. Note that B$ is 
zero along the rotation axis because the rotation speed 
goes to zero at 90° latitude. B$ is also zero everywhere 
along the equator because the radial component of the 
magnetic field goes to zero there, for a dipole, and it is 
the radial component of the magnetic field that is twisted 
to generate B$. Thus, for a twiste d dipole field, B ^ has a 
maximum value at mid latitudes ll Washimill 1 990(1 . There 
is a magnetic force associated with gradients in B^ (pro- 
portional to — V-B^), which point toward the rotation 
axis (as studied by, e.g.. ISakurailll985() and also toward 
the equator (as pointed out by WS93). For the ISO2000 
case, these magnetic forces are too small to have a sig- 
nificant effect, and the wind is accelerated hydrodynam- 
ically, though modified by the greater (less) than radial 
divergence near the polar (equatorial) region, imposed 
by the poloidal magnetic field. 

For the cases with more rapid rotation, however, the 
forces associated with B§ become important, and act 
both perpendicular to and parallel to the poloidal field 
to modify the flow. This is evident in the right panel of 
Figure ^ which illustrates the steady-state flow for the 
ISO100 case. The quantities shown in the right panel are 
the same as for the left panel, except that the greyscale 
image of B^, is scaled so that black corresponds to 860 
Gauss, and white is zero. The only difference between 
the two cases is that the star in the ISO 100 case rotates 
20 times faster than for ISO2000, so a stronger B^ is gen- 
erated that influences the wind. A comparison between 
the two panels of Figure ^ reveals this influence. First, 
the dead zone is significantly smaller for ISO100 because 
the magnetic pressure gradient in B^, acting perpendic- 
ular to the poloidal field, is strong enough to "squash" 
the closed field region. In addition, the field lines in the 
ISO100 case exhibit greater than radial divergence from 
the polar region of the star (as in the ISO2000 case) but 
then undergo a change in curvature outside a few i?* and 
begin to diverge more slowly than radially near the ro- 
tation axis. This gives rise to an asymptotic collimation 
of the flow, and is caused by the magnetic pressure gra- 
dients in B,f, that are perpendicular to the poloidal field 
and directed toward the axis. This collimation results 
in shallower pressure gradients along the pole, leading to 
slower acceleration of the wind there, and the sonic point 
is shifted outward relative to the ISO2000 case. This ef- 
fect (and due to the more slowly diverging poloidal mag- 
netic field) also moves the Alfven surface further out on 
the poles. The most important difference between the 
left and right panels of Figure Q — that a rotating dipole 
field simultaneously compresses the flow on the magnetic 
equator and collimates it along the rotation axis — was 
first shown by WS93 (e.g., compare with their fig. 1). 

Fast rotation also affects the flow in the direction par- 
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Fig. 2. — Forces parallel to the reference field line indicated by the 
light grey line in right panel of Figure as a function of position 
s along the field line. Lines labeled P, G, C, and M represent the 
pressure gradient,gravitational, centrifugal, and magnetic forces, 
respectively feqs. latSl . The locations s s and sa indicate the sonic 
and Alfvenic critical points, respectively. 



allel to the poloidal magnetic field. Figure [21 shows these 
forces along a reference magnetic field line that connects 
to the star at 65° latitude, for the ISO100 case. This 
field line is indicated by the light grey line in the right 
panel of Figure and the data in Figure [21 is from the 
second computational box. The forces per mass a long a 
field line are given by fe.g.. lUstvugova et al.111999]) : 



fp 

fo 
fc 

/m : 



1 dP 
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where the subscripts P, G, C, and M correspond to the 
pressure gradient, gravitational, centrifugal, and mag- 
netic forces, respectively. Here, s, is the distance along 
the poloidal field, G is the gravitational constant, and 
is the mass of the star. Also note that we have used 
both the cylindrical (r) and spherical [R) radial coordi- 
nate, where appropriate. Figure [21 shows these forces in 
units of cm s -2 (for r\ Car parameters, where the surface 
gravity is 270 cm s _1 ), as a function of position along s, 
where li?* corresponds to the stellar surface. It is evi- 
dent that the pressure gradient force along the reference 
field line dominates near the star, but the centrifugal 
and magnetic forces are significant. We find that the 
ratio of fc/ fv reaches a maximum value of ~ 30% at 
~ 10-R* along s. Consequently, material at mid latitudes 
is partially magnetocentrifugally accelerated, leading to 
a faster outflow. As a result, the sonic and Alfvenic sur- 
faces are closer to the star than for the ISO2000 case (see 
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Fig. 3. — Greyscale of (white is zero, black is strong) with selected magnetic field lines overdrawn in grey for two cases where the star 
has a high temperature polar cap. The light grey field line in each panel intersects the stellar surface at 74°, the latitude of the temperature 
discontinuity on the star. The dashed line indicates the radial Alfven surface, and the dotted line is the sonic surface. The left panel 
contains data for the very weakly rotating case (T = 2000 days), while the right panel corresponds to faster rotation (T = 100 days). The 
dash-dotted line in the left panel is the location of the sonic surface for a simulation without rotation, and without a magnetic field. 



Fig. Hi . 

The left and right panels of Figure [3] show the steady- 
state winds of the HTC2000 and HTC100, respectively, 
in the same format as Figure ^ (the grey-scaling for B$ 
is also the same for each respective panel). The two 
cases shown in Figure |3 are identical to those in Figure 
Q except that the star now has a high temperature polar 
region. In each panel of Figure El there is a light grey 
field line that is anchored on the stellar surface at 74° 
latitude. Since that is the latitude of the edge of the 
high temperature region, and since the field line is also a 
streamline, it marks the division between material that 
originates from the hot polar cap and material that flows 
from lower latitudes. 

A comparison of the left panels of Figures ^ and [3] 
reveals the effect of the hot polar cap on the wind from 
a slowly rotating star with a dipole field. The increase in 
thermal energy on the pole more quickly accelerates the 
wind, which moves the sonic and Alfven surface much 
closer to the star near the pole. However, since the hot 
polar wind is not in pressure equilibrium with the wind at 
lower latitudes, it expands rapidly toward the equator as 
it flows outward (more than for ISO2000), until there is a 
meridional balance in pressure between the high and low 
latitude flows. Far from the star, the flow in a large range 
of latitudes originates from the hot polar cap (further 
discussed in ij4 . 3fl . The convergence of flow toward the 
equator reduces radial pressure gradients there, which 
moves the sonic and Alfven surface outward (relative to 
the ISO2000 case). The azimuthal magnetic field exhibits 
the same qualitative behavior in both the ISO and HTC 
cases. 

A comparison between the left and right panels of Fig- 



ure |3 reveals the effect of significant rotation in the HTC 
models. As with the ISO cases, centrifugal flinging accel- 
erates the flow more rapidly in the equatorial direction, 
moving the location of the sonic and Alfven surface in- 
ward. Also, the dead zone is somewhat squashed by the 
magnetic forces (associated with B^) directed toward the 
equator. However, in contrast to the ISO cases, the dy- 
namical effects of the rotation in the HTC case are not 
as important at high latitudes, due to the stronger pres- 
sure gradient there. Along a field line connected to the 
HTC100 star at 80° latitude, for example, /c//p reaches 
a maximum value of less than 1%, compared to a ratio of 
greater than 9% along a field line connected to the star 
at 65°. This difference between the high and low lati- 
tude flow is evident in the right panel of Figure |3 as the 
lowest latitude open field line is bending slightly pole- 
ward (collimating), while the higher latitude field lines 
are not. The high latitude flow is is inherently more en- 
ergetic, so it would require a stronger field and/or faster 
rotation than the lower latitude flow, to be significantly 
altered. Consequently, the effect of significant rotation 
on the HTC models is to "clear out" the wind at mid 
latitudes as it is compressed toward the equator and also 
against the side of the polar flow. Note that the direc- 
tion of the light grey line points to higher latitudes in 
the HTC100 case than in the HTC2000 case, indicating 
that, when the star spins faster, the energetic polar flow 
is confined to a narrower opening angle. 

4.2. Global Outflow Properties 

The presence of a dipole field and the rotation of the 
magnetized star influences the global outflow rates of 
mass (M w ), angular momentum (J w ), and energy (E w ). 



10 



Matt & Balick 



£2 = 



T = 200 days 



T = 100 days 



T = 50 days 




O -3 






90 60 30 90 60 30 90 60 30 90 60 30 



600 



400 



200 




600 



400 



200 




600 



400 



200 ; 




600 



400 



200 




90 60 30 
9 (degrees) 



90 60 30 
9 (degrees) 



90 60 30 
9 (degrees) 



90 60 30 
9 (degrees) 



Fig. 4. — The mass outflow rate, 8M /<9(sin$), (top panels) and radial velocity of the wind (bottom panels) as a function of latitude, 9, 
for cases with (from left to right) no rotation, T = 200, T = 100, and T = 50 days. The dashed lines represent cases with spherical thermal 
wind driving, while the solid lines are for the cases with a high temperature polar wind. All panels are calculated at a distance of 40-R* 
from the star. 



We calculate these quantities from the simulation data 
by integrating the fluxes over an enclosed spherical sur- 
face at a radius R. Thus, 



M w = AttR 1 



J w = AttR 2 



pvn d(sin0), 



PVr A d(sm 9), 



(6) 



(7) 



E w = AttR 2 / pv R E' d{sm 9) , 



(8) 

where 9 is the latitude, vr is the velocity component in 
the spherical radial direction (vr = v ■ R), and 



A = v^r 



E' 



7- 



7 P _ GM* 
1 p R 
VtpB^Bp 



(9) 



(10) 
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Here, A is the integral of motion discussed in section ET4l 
that gives the angular momentum per mass carried in 
the wind, and E 1 = E + Aft gives the specific energy 
carried by the wind. The first three terms in equation 
E3represent the kinetic (£k), thermal, and gravitational 
energy, while the last two terms together represent the 
magnetic energy (Em = Poynting flux divided by pvr). 
Note that the equations have been multiplied by a 
factor of two to take both hemispheres into account. 



The top row of panels of Figure 0] shows the differ- 
ential mass loss rate, dM w / 'd(sm8) = 2nR 2 pvfi ("mass 
flux"). The bottom row of panels shows the radial veloc- 
ity vr. Both the mass flux and vr are plotted versus 9, 
at a constant spherical radius of 40i?» (which equals 300 
gridpoints in the third computational box), for several 
models. From left to right, the panels for both mass flux 
and vr show cases with increasing rotation rate. In each 
panel, the dashed line is for the ISO case, while the solid 
line represents the HTC case. Note that for the ISO2000 
and HTC2000 cases, which are featured in the left panels 
of Figures H an d El the rotation is so slow that the mass 
flux and vr for those models are the same (within a few 
percent) as the data plotted in the left panels of Figure^] 
for the ISOoo and HTCoo cases. Effectively, then, Figure 
0| displays data from all 10 of the models with non-zero 
magnetic field. 

The dashed lines in the left panels of Figure 01 reveal 
that steady-state wind in the non-rotating ISOoo case is 
not completely isotropic. For reasons discussed in sec- 
tion 14.11 th e flow is densest near the equator (see also 
iMatt et alJ 1200(1 and AO} . This is the effect of the 
dipole field on an otherwise isotropic wind. For the HTC 
models, the flow is generally faster at higher latitudes, 
and the mass flux is higher than for the ISO cases, due 
to the more energetic wind driving in the HTC cases. 

The integrated outflow rates M w , J w , and E w are listed 
in the second, third, and fourth columns of tabic for 
each case. Since E w is dominated by the thermal en- 
ergy within the computational domain, due to our chosen 



Single- Wind MHD Model for rj Car 



11 



TABLE 2 
Simulation Results 
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a Given is the value at a radius of 80/?*. 
b Given is the value at a radius of 40i?*. 



value of 7 near unity, we list the individual components 
of Ek and Em in the fifth and sixth columns of table 
121 These are calculated with equation [H] but using only 
the kinetic or magnetic terms of equation 1101 The table 
lists these values at a spherical radius of 80i?* (which 
equals 300 gridpoints in the fourth computational box), 
where the gravitational energy is negligible. If one scales 
M w and i?* down to solar values for the ISO2000 case 
(as described in M3.3[) . this case predicts J w « 1.1 x 10 30 
erg for the solar wind. This value is comparable to the 
predictions of J w ~ 10 30 erg from the solar models of 
iKeppens fc Goedbloedl 1)2000(1 . and this further verifies 
our numerical solution. 

It is evident from tabled that M w and _E W are slightly 
less in the ISOoo case than for ISO°oo (and similarly for 
HTCoo and HTC°oo). This is because the existence of 
the dead zone "shuts off" part of the flow. However M w 
and E w are not simply reduced by the amount propor- 
tional to the fractional area covered by the dead zone on 
the star, since a reduced wind flux near the equator is 
counteracted, somewhat, by an increased flow from the 
poles, driven by the faster than radial divergence of the 
wind there. In order to quantify the anisotropy of the 
wind, we use a crude measure given by the ratio of the 
maximum value of vr to the minimum vr found in the 
wind at a spherical radius of 40i?*. This ratio is listed 
in the seventh column of table for all cases, while the 
eighth column contains a similar ratio for the anisotropy 
in the mass density. 

From the data in table El and plotted in Figure 01 it is 
evident that, as the stellar rotation rate increases, M w , 
J w , E w , and the anisotropy of the flow increases. The 
effect of strong azimuthal magnetic field is to clear out 
material from mid latitudes, where it has a maximum 
strength, as it both compresses material toward the equa- 
tor, forming an outflowing disk, and collimates material 
toward the pole. For the ISO models, a jet appears as an 
increase in mass flux along the rotation axis, as evident 
in the top row Figure^ For the HTC models, the open- 
ing angle of the energetic flow from the pole of the star 
is decreased an outflowing "lobe" appears — as rotation 
increases. These features are discussed further in section 



PI Note that for the ISO50 and HTC50 cases, the mass 
fluxes and vr are similar at latitudes less than roughly 
40°. 

Also evident in Figure 01 is that, in each case, there is 
a general anti-correlation between the mass flux and the 
wind velocity. That is, where the mass flux is high rel- 
ative to the non-rotating case, the velocity is relatively 
low. This is because the convergent flow that increases 
the local mass flux also makes the thermal pressure gra- 
dients shallower, leading to a slower acceleration of the 
flow. However, in spite of the general decrease in vr near 
the rotation axis and equator, we find that the linear mo- 
mentum and kinetic energy flux is generally the largest 
where the mass flux is largest. 

An exception to the mass flux and wind velocity anti- 
correlation occurs for the ISO50 case (dashed line in the 
upper right panel of Fig. 0J), where the mass flux is en- 
hanced near the rotation axis (as in other cases), but 
decrease to a minimum value on the axis, resulting in a 
"hollow" jet structure. This behavior is similar to cold, 
magnetocentrifugally driven winds, i n which the inner 
part of the jet is generally hollow fe.g.. lOuved fc Pudritzl 
Il997|) . The flux in the center of this jet is small be- 
cause the magnetic field cannot centrifugally fling ma- 
terial along the exact rotation axis, where the rotation 
speed goes to zero. In our models, we find that, along 
a field line that connects at 9 = 65° on the star, /c//p 
reaches a maximum value of 10%, 30%, and 70% for the 
ISO200, ISO100, and ISO50 cases, respectively. So the 
sequence of ISO models is showing a transition from a 
purely thermal wind, to a thermo-centrifugal wind (see 
WS93). This is also evident in table |3 where the ratio 
of Em/E w increases with the rotation rate. 

4.3. Morphology: Disks, Jets, and Bipolar Lobes 

The outermost (fourth) simulation box reaches beyond 
100i?», and the conditions in the steady-state wind there 
are indicative of the shape of the flow very much further 
from the star. The three panels of Figure show loga- 
rithmic density contours from the fourth box for three 
ISO cases with different stellar rotation rates (from left 
to right: ISOoo, ISO200, ISO100). For reference, the 
thickest contour line in each panel corresponds to a den- 
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r (R.) r (R») r (R») 

Fig. 5. — Logarithmic density contours in the r-z plane for cases with (from left to right) no rotation, T = 200, and T = 100 days and 
where the star has a constant temperature across its surface. Shown is data from the outermost computational grid, and the star is at the 
origin. The thickest contour line corresponds to a density of 1.7 X 10~ 13 g cm -3 , and the spacing between each contour is a factor of two. 
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Fig. 6. — Same as Figure|S] but for cases with (from left to right) no rotation, T = 100, and T = 50 days and where the star has a high 
temperature polar cap. The dashed line in each panel shows the magnetic field line that intersects the stellar surface at 74° , the latitude 
of the temperature discontinuity on the star. 



sity of 1.7 x 10 -13 g cm~ 3 . The left panel shows the 
effect of the dipole field alone — the wind is denser near 
the equator than along the pole. For stars that rotate 
quickly (middle and right panel), notice the formation of 
the jet along the rotation (z) axis and the disk along the 
equator. Also, compare the middle panel with figure 2 
of WS93. 

Figure is the same as Figure |S] but for the HTCoo, 
HTC100, and HTC50 models (from left to right). The 
left panel shows the conditions in the wind from a star 
with both a dipole magnetic field and an energetic po- 
lar wind. The density is enhanced along the equator, 
due to the dipole field, but also along the poles from 
the increased thermal driving there. The dashed line in 
each panel follows the magnetic field line that originates 
at a latitude of 74° on the stellar surface (as the light 
grey lines in Fig. So the material that exists at a 
higher latitude than the dashed line originates from the 
high temperature region on the star. In the left panel 



(where the star does not rotate), only material at very 
low latitudes originates from the cooler region on the 
star. When stellar rotation becomes important (mid- 
dle and right panels), the mid latitudes are cleared out 
as magnetic pressure gradients compress material toward 
the equator and against "walls" of the polar flow. Instead 
of a narrow jet, wide-angle, bipolar lobes are formed. Ul- 
timately, the opening angle of the lobes is determined by 
the balance between the thermal pressure in polar flow 
and magnetic pressure (associated with in the low- 
latitude flow. 

For the HTC50 case (right panel of Fig. [fjj), the den- 
sity contours superficially resemble 7/ Car's homunculus. 
The morphology of the flow is more clearly seen in the 
projected, 3-dimensional, isodensity surfaces plotted in 
the top row of Figure On the top left is an isoden- 
sity surface in the steady-state wind of the ISO50 case, 
and the top right is from the HTC50 case. The surface 
on the top right corresponds to the same density as the 
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Fig. 7. — Isodensity surface in the steady-state wind for two 
different models (top), each with a stellar rotation period of 50 
days, compared to two HST/WFPC2 images of nebulae (bottom). 
Top left: density surface at 5.2 X 10 — 14 g cm -3 for a case with 
isotropic thermal wind driving. Top right: density surface at 1.7 X 
10 -13 g cm -3 for a case with an enhanced polar wind. These are 
projections with the rotational symmetry axis tilted 30° from the 
plane of the sky. Bottom left: proto-planetary nebula IRAS 17106- 
3046 lKwok "etaT]l20001) . Bottom right: r\ Carinae IMorse et alJ 
1998). 



thick line in the right panel of Figure which serves as 
a size reference for the 3-D projections of Figure [7| To 
generate these images, we have exploited the symmetries 
in the system by rotating the 2-D density contour from 
the simulation about the rotation axis and reflecting it 
across the equator. Notice that the disk is similar for the 
two cases, but the polar flows are drastically different. 
For comparison, we have also included HST images of 
two nebulae on the bottom row of the Figure. The bot- 
tom l eft is the proto-planetary neb ulae IRAS 17106-3046 
(from lKwok. Hriv nak, fc Sul|2000t), and the bottom right 
is 77 Car (from IMorse et all 119981) . Both nebulae have 
their symmetry axes along a line from bottom left to top 
right, both exhibit equatorial structures, and both are 
oriented such that the axis pointing toward bottom left 
of the image is tilted out of the plane by (very roughly) 
30°. 

The top right panel of Figure bears a remarkable 
resemblance to the skirt and h omunculus morpholog y 
of rj Car (bottom right; also see lDavidson et all feOOl). 
However, the comparison is only superficial because we 
have considered only the conditions in a steady-state 
wind, while the lobes and skirt of rj Car were ejected 
in an outburst. Also , the lobes of rj Car are hollow 
ijDavidson et alJl2001j) . It is not clear from this prelim- 



inary study what sorts of morphologies would result, if 
the wind were, for example, "turned on" for a few years 
(simulating an outburst) and then "turned off" again af- 
terwards and left to drift for 150 years. Such a study is 
left for future work. 

5. SUMMARY AND DISCUSSION 

We have studied a thermally-driven wind from a star 
with a rotation-axis-aligned dipole magnetic field and 
varying amounts of rotation. We have considered both 
isotropic thermal driving and that which is enhanced on 
the pole of the star. Using MHD simulations, we deter- 
mined the 2.5-dimensional (axisymmetric), steady-state 
wind solution. Our parameter study resulted in winds 
with a wide range of shapes. Winds are slightly enhanced 
along the magnetic equator by the presence of a dipole 
magnetic field. When the star rotates, the poloidal mag- 
netic field is twisted azimuthally, and magnetic pressure 
gradients associated with direct material away from 
mid latitudes — toward the rotation axis and toward the 
magnetic equator. For cases with isotropic thermal driv- 
ing, fast stellar rotation results in an outflowing disk and 
jet morphology. For cases with an enhanced polar wind, 
there is an outflowing disk and wide-angle bipolar lobes. 
As far as we know, no physical model for producing such 
a wind (with a simultaneous disk and lobe morphology) 
has yet been demonstrated in the literature. 

5.1. Implications for rj Car 

The disk and lobe morphology of our steady-state wind 
solutions resembles the skirt and homunculus morphol- 
ogy of the 77 Car nebula, though the latter was formed 
in a single eruptive outburst and is much more complex 
than the flows of our study. The opening angle of the 
lobes in our model HTC50 most resembles that of the 
homunculus. In addition, the total energy outflow rate 
in that model is ~ 9 x 10 40 erg s _1 (table EJ), which is 
comparable to the observed lu minosity of 77 Car during 
outbu rst of ~ 8 x 10 40 erg s _1 fDavidson & Hum phreys! 
1997). For a comparison between the HTC50 model and 
the 77 Car nebula to be appropriate, 77 Car needs to have 
had a roughly axis-aligned dipole magnetic field of the 
order of 2.5 x 10 4 Gauss during the great eruption in the 
mid 1800's. It also needed to have a wind driving mech- 
anism that produced a flow that was inherently more 
energetic on the poles and needed to be rotating at a sig- 
nificant fract ion of breakup speed (both o f which are also 
predicted by Dwarkadas & Owocki 2002). Assuming, for 
now, that 77 Car's wind resembled the HTC50 model dur- 
ing the great eruption, our model implies several things. 

First, we find that ss 30 G at 40i?* and at mid lat- 
itudes. Assuming r~ 4 dilution of the field and a size of 
37500 AU for the 77 Car nebula ijDavidson &: Humphreys! 
Il997|) . this predicts a ~ 15 mG (azimuthally directed) 
field in the present-day ho munculus. This is la rger than 
the ~mG fields detected by Ai tken et alJ l|1995[K but con- 
sidering the uncertainties in our assumed parameters, 
and that the magnetic field can be dissipated by other 
processes, this comparison should not be stringent. An- 
other unique prediction of magnetic wind theory is that 
the entire flow should be rotatin g (as verified for proto- 
stellar iets: lBacciotti et al.l2002D . The HTC50 model has 
a rotational speed of ~ 30 km s _1 at 40i?* in the disk. 
Assuming angular momentum conservation in the wind 
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from that point outward, this predicts a rotation speed 
of — 20 m s _1 in the skirt at 25000 AU, which would be 
difficult to detect observationally. 

The rotation of the outflow extracts angular momen- 
tum from the star, and we find a loss rate of ~ 1.3 x 10 46 
erg (tableEJ during outburst. The angular momentum of 
r\ Car, for the assumed parameters, is J* = k 2 M^Sl^R 2 m 
1.4 x 10 55 fc 2 erg s (where k is the normalized radius of 
gyration). Assuming k 2 ~ 0.1-0.2, the total amount of 
angular momentum extracted by the wind during the 
outburst would be comparable to J* in only 3.5-7 years! 
The reason for such a large J w is mainly because of the 
large mass of the homunculus, which is > 1% of the mass 
of the star. Particularly, if this much mass is shaped mag- 
netocentrifugally, as we have considered here, J w will 
necessary be large. This result suggests either that rj 
Car can have only one or two such outbursts during its 
lifetime, or that the star must gain angular momentum 
from an external source. No te that 77 Car may have a 
binary companion (see, e.g.. ISmith et al.ll2004f) . which 
could likely further influenc e the structure of the wind 
and possibly spin up 77 Car ijSmith et alJl2003|) . though 
a detailed calculation remains to be done. 

5.2. General Application 

The range of morphologies exhibited in our parame- 
ter study may be applicable (with appropriate modifica- 
tions) to a wide range of observed outflow nebulae. A 
general conclusion of our study is that, for any magnetic 
rotator wind in which collimation occurs, if the source 
has an aligned dipolar field, the formation of an outflow- 
ing disk is inevitable. So any outflow nebulae that ex- 
hibits quadrupolar symmetry (i.e., outflowing disk plus 



jets or disk plus lobes) may be explained with winds 
similar to those presented here. For example, the gen- 
eral torus pl us jet morphology of the crab pulsar x-ray 
nebula (e.g.. lHester et ai1l2002j) could be explained by a 
wind accelerated in a rotating dipole field, though rel- 
ativistic effects and other pulsar physics should be con- 
sidered. In addition, proto-planetary nebulae all show 
axisymmetries, and the overwhelming majority of opti- 
cally visible proto-planetary nebulae have disks seen in 
silhouette (e.g., see IRAS 17106-3046 in the lower left 
panel of Figured). A mechanism that generates disks 
and lobes together may turn out to be important. The 
general model presented here is also appealing because of 
the relatively small number of parameters, and its ability 
to link observed parameters to the conditions very near 
the source of the wind. 

It is important to note, however, that many astro- 
physical nebulae of various types are evidently shaped 
by processes that are inherently time-dependent, either 
through the interaction of multiple winds or in outburst- 
type events. Time-dependence adds many additional pa- 
rameters to any outflow model, making parameter stud- 
ies challenging, but greatly enriching the range of pos- 
sible outcomes. Future work with magnetocentrifugal 
stellar winds should include time-dependent effects. 
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